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Abstract 

In this paper, we apply results on number systems based on continued fraction 
expansions to modular arithmetic. We provide two new algorithms in order 
to compute modular multiplication and modular division. The presented algo- 
rithms are based on the Euclidean algorithm and are of quadratic complexity. 



1. Introduction 

Continued fractions are commonly used to provide best rational approxima- 
tions of an irrational number. This sequence of best rational approximations 
{Pi/<li)i&i is called the convergents' sequence. In the beginning of the 20 th cen- 
tury, Ostrowski introduced number systems derived from the continued fraction 
expansion of any irrational a Q. He proved that the sequence (qi)i<=n of the 
denominators of the convergents of any irrational a forms a number scale, and 
any integer can be uniquely written in this basis. In the same way, the sequence 
(qia — Pi)i£ti also forms a number scale. 

In this paper, we show how such number systems based on continued fraction 
expansions can be used to perform modular arithmetic, and more particularly 
modular multiplication and modular division. The presented algorithms are 
of quadratic complexity like many of the existing implemented algorithms 0, 
Chap. 2.4]. Furthermore, they present the advantage of being only based on 
the extended Euclidean algorithm, and to integrate the reduction step. 

In the following, we will first introduce notations and some properties of the 
number systems based on continued fraction expansions in Section [5] Then we 
describe the new algorithms in Section[3] Finally, we give elements of complexity 
analysis of these algorithms in Section |4l and perspectives in Section 



2. Number systems and continued fractions 

2.1. Notations 

First, we give some notations on the continued fraction expansion of an 
irrational a with < a < 1 [3]. We call the tails of the continued fraction 



expansion of a the real sequence (rj)jgN defined by 



r = a, 

n = 1/rj-i - [l/fV-iJ. 

We denote (fci)igN the integer sequence of the partial quotients of the continued 
fraction expansion of a. They are computed as fej = [l/^-iJ- We have 

a = := [0; fci, fc 2 , • • • , h 



hi 



hi + rj 

We write Villi the i th convergent of a. The sequences (pi)igN and (gi)igN 
are integer valued and positive, 

— = [0;fci,fc 2 ,...,fcj]. 

We will also write (6>i) i6 N the positive real sequence of (— l)*(<fca — pi) which 
we call the sequence of the partial remainders as they are related to the tails 
by Ti = 9i/6i-i. Hereafter, we recall the recurrence relations to compute these 
sequences, 

P-! = 1 po = Pi = Pi-2 + kiPi-1, 

q-i = q a = 1 qi = <? 4 _ 2 + fc,gi_i, 
0_i = 1 6>o = a 0j = #i_2 — kiOi-%. 

We also write 77, = g^a — pi the sequence of the signed partial remainders, 
which elements are of sign (—1)*. The sequence (rji)^ of the signed partial 
remainders can be computed as ((— l) l 0i)i e N- 

2.2. Related number systems over irrational numbers 

In this section, we present two number systems based on the sequences of 
the signed partial remainders (ry^igN & n d the denominators of the convergents 
(<Zi)ieN of an irrational a. They have been extensively studied during the second 
part of the 20" 1 century 

Property 2.1 ([1, Proposition 1]). Given (qi)i<=n the denominators of the con- 
vergents of any irrational < a < 1, every positive integer N can be uniquely 
written as 

m 

N = 1 + y^n^_i 

i=l 

where 1 ^ — n * ~ ^ 1 ^'^ — n% ~ ^' * — ^' ("Markovian" conditions). 
I »i = 1/ = fc l+ i 
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Algorithm 1: Integer decomposition in Ostrowski number system, 
input : JVeN, (qi)i <m 

m 

output: rii such that N = 1 + ^] n%qi-\ 

i=i 

1 imp A^ — 1 ; 

2 i(-m; 

3 while i > 1 do 



m <- [tmp/qi-i\; 
tmp <— tmp — riiqi-i] 
i <— i — 1: 



This number system associated to the (<2i)ieN is named the Ostrowski number 
system. To write an integer in this number system, we use a classical decompo- 
sition algorithm (Algorithm [TJ . The rank m is chosen such that q m > N. 

Property 2.2 ([1, Proposition 2]). Given (jji)ieN the sequence of the signed 
partial remainders of any irrational < a < 1, every real ft, with < f3 < 1 
can fee uniquely written as 

/8 = a + bifji-i 

i=l 

where ( ? ^ * ^ ** " L J <h< h, for i > 2, 

There also exists two other number systems that are dual to these two. One 
decomposes integers in the basis ((— l) l qi)i^n and the other decomposes reals in 
the basis of the unsigned partial remainders Q. The second Markovian 

condition then becomes 6j+i = if hi = fcj. An algorithm to write real numbers 
in the (0j)jerj number scale has been proposed by Ito @. It proceeds by iterating 
the mapping T x : (a, (3) -4 (1/a - |l/a| , /3/a - L/V"J )• 

2.3. Related number systems over rational numbers 

In this subsection, we consider a = p/q rational. We recall that the continued 
fraction expansion of a rational is finite. We denote 

P 

- = [0;fci,fc 2 ,---,fcn] 

the continued fraction expansion of p/q, and recall p n — p and q n = q. 

The Ostrowski number system still holds for integers N < q ni since the 
keypoint in the Ostrowski number system is that there exists q rn such that 
q m > N. 

The {rji)i <n number system also still holds under one supplemental condition: 
j3 must be rational with precision at most q (i.e. the denominator of (3 must be 
less or equal than q). 
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3. Modular arithmetic and continued fraction 



In this section, we consider a — a/d. We highlight that the same decompo- 
sition (61, ... , 6 n +i) can be interpreted in two ways depending on the number 
system used. In the Ostrowski number system, we obtain an integer N whereas 
in the number scale (?7i)igNj we obtain the reduced value of Not mod 1 
Hence, we will use the fact that studying an integer a modulo d is similar to 
considering the rational a/d modulo 1. This enables us to use properties 12.11 
and 12.21 to compute modular multiplication and division. 

3.1. Modular arithmetic and continued fraction 

First, we briefly recall how continued fraction expansion and the Euclidean 
algorithm are linked. We write (^)ign the integer sequence of remainders when 
computing gcd(a, d). This sequence is composed of decreasing values less than 
d. We also write (f^)igN the sequence ((— l)^-)^- We obtain the following 
recurrence relation, and recall the recurrence relation over the (f^gpj sequence 
of partial remainders of the continued fraction expansion of a/d : 

Q'-i = d O'o = a 6'i = O'i-2 - Wi-i/Oi-iM-i 
0-i = 1 #0 = a/d 9i = #i_2 — \ 9i-i/9i-i\9i-i- 

It is widely known and can be easily proved by induction that both sequences 
compute the same partial quotients, that we will note fcj. 

Proof of k l+1 = L^-iMJ = L^-i/^J- Weproveitbyprovmg^_ 1 /e i = ^_ 1 /^. 

• Base case : Q-\/Qq = d/a = 0'_ 1 /6' 

• Induction : Let i such that Qi-i/Qi = O'i-i/Q'i- 

Oi-i 6'i-i 



which implies Bi/6 i+1 = 0'J6' i+l . □ 

It can also be noticed that 77- = rjid. Actuality, B' i = 6 id as the extended 
Euclidean algorithm compute the relations d\ = (—\) l (qia — pid). In particular, 
it gives the Bezout's identity with 9' n _ 1 = (— l) n_1 (g„_ia — p n -id) = gcd(a,d), 
and q n -i the inverse of a if a is invertible modulo d (gcd(a, d) = 1). 
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3.2. Modular multiplication 

Now, given a,b G Z/eK, we write c = a-b mod d the integer < c < d such 
that ab — [ab/d\ ■ d = c. 

We can observe that the decompositions presented in properties 12.11 and 
12.21 are both unique and both need the same "Markovian" condition over their 
coefficients. Hence, we can interpret the same decomposition in both basis. 

Theorem 3.1. Given a, b G Z/dZ, and {qi)i< n , {Vi)i<n from Euclidean algo- 
rithm on a and d, if we write b in the (qi)i< n number scale as 

71+1 

b = 1 + ^ Mi-l> 

then 

n+l 

a ■ b mod d = a + biii' i _ 1 . 

i=l 

Proof. First, we consider b < q n , it can be written in the Ostrowski number 
system as 

n 

b= i + y^bjgi-i, 

i=l 

and the coefficients bi respect the "Markovian" condition of the Ostrowski num- 
ber system. Hence, 

n 

a ■ b = a + biqi^ia. 

i=l 

By definition, r\i — qia — pi, thus 

n n 

a-b = a + ^ biVi-l + X] b iPi-l- 

i=l i=l 

As the coefficients b^s verify the "Markovian" condition, the uniqueness of the 
decomposition in property 12 . 2 1 gives < a + J27=i biVi-i < 1 and X)i=i biPi-i G 
N. Hence, 

n 

a ■ b mod 1 = a + frj?7i-i • 

i=l 

By multiplying this inequality by d, as a = a/d and rj^ = r/id, we obtain 
a ■ b mod <i = a + b; t ri' i _ l . 

i=l 

which finalizes the proof of the theorem for b < q n . 

Now if b > q n and b = b n+ \q n + b' with b' < q n the remainder of the 
division of b by q n , b' can be uniquely written in the Ostrowski number system. 
Furthermore, as rj' n — 0, b n+ irj' n — 0, which finishes the proof. □ 
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3.3. Modular division 

Inversely, given a, b E Z/dZ, with a invertible modulo d (gcd(a, d) = 1) we 
can efficiently compute a^ 1 ■ b mod d. 

Theorem 3.2. Given a, b G Z/dZ with gcd(a, d) = 1, and {qi)i< n , [&i)i<n from 
Euclidean algorithm on a and d. if we write b in the (^)j<« number scale as 

n+1 
i=l 

n+1 

i/ien i/ we denote c = l) 1 ^ 1 , 

i=i 

a -1 • 6 mod d 6 {c, d + c}. 

Proof. The proof of correctness is similar to the one of theorem 13.11 using the 
facts that = 6^d and that 0j = (— l) l (%a< — _p;)- 

Now, the greatest integer c is clearly the one associated to the decomposition 
(ki, 0, /c3, 0, . . . , k n ) when n is odd. However, kiqi-x — qi — qi- 2 by definition, 
which implies 

(n-l)/2 

^ k 2i+1 q 2i = q n - 
i=0 

The smallest integer that can be returned is clearly the one associated to 
the decomposition (0, k 2 , 0, k^, . . . , k n ) when n is even. Once again, as fc^-i = 
q% ~ qi-2, we get 

n/2 

-^k 2l q2i-i = 1 - q n - 
i=i 

Hence, — d < J27=i l) t ~ 1 qi-i < d, that is to say, the result needs at 
most a correction by an addition by d. □ 

We mention that we also tried to decompose b in the (r/'^iKn signed re- 
mainders number scale and evaluate this same decomposition in the {qi)i< n 
number scale to compute modular division. We used Ito T 2 transform [5] 
T 2 : (a,f3) — s> (1/a — LV a Jj l7V a l ~ P/ a )- in practice, it returns the right 
result without the need of any correction. However, as the decomposition com- 
puted by Ito T 2 transform does not verify the same "Markovian" conditions as 
in the Ostrowski number system, we were not able to give a theoretical proof 
that it always returns the reduced result of the modular division. 

4. Elements of Complexity Analysis 

In this section, we introduce elements of complexity analysis of the proposed 
modular multiplication algorithm based on theorem 13.11 The same analysis 
holds for the division. 



G 
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Max expected 6„+i 

Fi gure 1: Probability law of the value of the coefficient 6 n _|_i 



First, the algorithm computes {qi)i< n and (rji)i<n- This can be computed us- 
ing the classical extended Euclidean algorithm in 0(log (d) 2 ) binary operations. 
We notice here that the divisions computed in the Euclidean algorithm can be 
computed by subtraction as the mean computed quotient equals to Khinchin's 
constant (approximately 2.69) [3j, p. 93]. Furthermore, big quotients are very 
unlikely to occur as the quotients of any continued fraction follow the Gauss- 
Kuzmin distribution 0, p. 83] 0, p. 352], 



P(fci = *) = -log 3 1 



(k + iy 



Second, the decomposition in (qi)i< n as in algorithm Q] also clearly has com- 
plexity in 0(log (d) 2 ). By the same arguments, the coefficients of the decompo- 
sition in (qi)i<„ can be computed by subtraction as they are likely small. The 
only quotient not following the Gauss-Kuzmin distribution is the coefficient & n +i 
as it corresponds to the quotient [b/q n \. We prove in Appendix A| that if a, d 
are uniformly chosen integers in [l,iV] and b is uniformly chosen in then 
when N tends to infinity, P(b n +i < k) tends to 



C(2)" X 



fc+i 

E 

,i=i 



(k + 1) 



(fc + l)C(3) 



Figure Q] shows the probability distribution of V(b n +i < k). In particular, 
we obtain P(6„+i < 3) w 92.5%. 

To finish the complexity analysis, evaluating the sum to return the final 
result can also be done in 0(log (d) 2 ). 
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5. Perspectives 

In this paper, we presented an algorithm for modular multiplication and 
an algorithm for modular division. Both are based on the extended Euclidean 
algorithm and are of quadratic complexity in the size of the modulus. 

Furthermore, the two stated theorems imply that, knowing the remainders 
generated when computing the gcd of a number a and the modulus d, one can 
compute efficiently reduced multiplications by a or a -1 . This can be useful 
in algorithms computing several multiplications and/or divisions by the same 
number a, as in the Gaussian elimination algorithm for example. 

The presented algorithms can also be useful in hardware implementation 
of modular arithmetic. They allow to perform inversion, multiplication and 
division with the same circuit. 

Further investigations have to be led to find optimal decomposition algo- 
rithms, that minimize the number of coefficients of the produced decomposition 
and their size. Also, we are working on an efficient software implementation of 
these algorithms. 
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Appendix A. Detailed proof of the distribution function of < k}. 

Let Ui, Ui and U3 be three independent uniform distributions over [0, 1]. We 
write a = \UiN], d = \U 2 N] and b = \U 3 d\. We denote A = {b < (k + l)q n }, 
B = {gcd(a, d) < k + 1}, B = {gcd(a, d) > k + 1} and B t = {gcd(a,d) = i}. 
Hence using the law of total probability we have 

F(A) = F(A DB)+ F(A n B), 

= |J F(ADB l )+ [_\ F(A n Bi), 

i<k+l i>k+l 

= \_\¥(A\B i )-¥(B i )+ \_\¥(A\Bi)-F(Bi). 

i<k+l i>k+l 

As the Bi are disjoint events, we have 

fe+l +oo 

F(A) = V(A\Bi) ■ P(Bi) + • P(fli). 

i=l i=k+2 

First, ¥(A\Bi) = lfori<fc + las6<d = gcd(a, d) ■ q n < {k + 1) • q n . 
Hence, 

fe+l +oo 
i=l i=k+2 

Now we want to determine F(A\Bi) for i > k + 2. Hereafter, we write 
Qi(-) =F(-\Bi) and 

F(A\B i ) = Q i (A), 

N N 

= 2 S = n {d = m}) ■ Q 4 (A I {a = ?} n {d = to}). 

1=1 m=l 

However, 

Qi(A | {a = Z}n{d = TO}) = — 

as 6 is uniformly distributed between 1 and d = iq n . If we consider the segment 
of length d and slice it in i segments of length q n , it can be interpreted as the 
probability that b is in the first k + 1 slices. Hence 

N N 

¥{A\Bi) -EE = n i d = m » ' 

Z=l m=l 

;=i m =l 
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As {a — 1} and {d — m} are independent by hypothesis (Ui and Ui are 
independent), 



,({a = /} n {d = m}) = Q,({a = /}) • = m}), 



and 



P(A|Bi 



fc + 1 



N 



N 



^Q l ({a = /})-^Q l ({d = m}). 



Z=l m=l 

Now, we use the fact that the sum of the probabilities over the whole sample 
space always sum to 1 to obtain 



K A \ B i) 



k + l 



If we recapitulate, 



fe+l 



fc+i 



i=l 



i=k+2 



Finally, it is widely known that P(-Bj) tends to — p 
infinity 0, p. 353]. Hence, we get 



C(2)~ 1 

— when N tends to 



fe+l 



lim P(A) = ^ 



cmr 1 



N-t+oo 



E 

i=k+2 



k + l C(2) _1 



C(2) _1 



'fe+l -. +oo 

Es + (*+d E * 



i=fe+2 



which equals to 



cm- 1 



'k+l /+co 1 fe+l 



= C(2)" 



fe+i . 



E 

.»=i 



i-(k + l) 



i 

+oo 



Ei 



\i=l 



By definition, Riemann zeta function equals 

+oo 



c(*) = E^ 



Hence we get the following simplification, which is more convenient for compu- 
tation and has been used to generate Fig. [IJ 



hmP( J 4) = C(2)- 1 

JV— >+oo 



fe+l 



,i=l 
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